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We show that for pulse coupled oscillators a class of phase response curves with both excitation 
and inhibition exhibit robust convergence to synchrony on arbitrary aperiodic connected graphs 
with delays. We describe the basins of convergence and give explicit bounds on the convergence 
times. These results provide new and more robust methods for synchronization of sensor nets and 
also have biological implications. 
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Synchronization in systems of pulse coupled oscillators 
(PCOs) is a fundamental issue in physics, biology and 
engineering. Examples from nature include synchroniza- 
tion of fireflies J^, Josephson junctions neurons in 
the brain [3] and the sinoatrial node of the heart [3]. In 
addition to the general study of pulse coupled oscillators 
^ [S] , there has been recent interest in their use for syn- 
chronization in sensor-networks [6l [7]. However, many 
of the existing models of PCOs, in particular those for 
which one can prove analytical results, are limited by 
strong assumptions. In this paper we analyze an inter- 
esting class of PCOs for which one can prove robust con- 
vergence results on arbitrary aperiodic connected graphs, 
even with propagation delays and a non constant graph 
topology (such as when spatially embedded nodes are 
mobile). This class of PCOs is of particular biological 
relevance because it explicitly includes both inhibition 
and excitation in the phase response curve (PRC), much 
like the type II PRCs seen in nature [SHSl E 1^. Addi- 
tionally, these PCOs provide guidance for the design of 
engineered systems of PCOs; improving on the current 
technology, by providing theoretical bounds for robust 
convergence under propagation delays and covering more 
diverse topologies. 

Our analysis was motivated by our prior work which 
used machine learning and genetic algorithms to engineer 
PRCs which would converge under propagation delays. 
In that work |10j we found that such algorithms typi- 
cally generate a very particular variety of type II PRCs. 
As we show below, for engineering applications such as 
sensor net synchronization [B], these PRCs are superior 
to those typically used and allow for a precise analysis 
which appear to differ from the majority of the literature. 
Namely our analysis does not rely on linear stability, in- 
stead our convergence argument makes use of values of 
the PRC over the entire domain as opposed to derivatives 
of the PRC at a single point. The analysis also shows 
that precise normalization of inputs is not required to 
achieve synchronization with propagation delays, unlike 
that suggested by the analysis in [11) . 

To begin, we describe the general structure of a PCO 
model on an arbitrary directed graph under delays. 



There are n oscillators where oscillator i's state is de- 
scribed by 0i(t) G [0,1]. (pi evolves with natural fre- 
quency 0=1 and emits a pulse as its phase is reset from 
1 to 0. The pulse is received time r < .5 later by all the 
successors of i, S{i) (predecessors denoted P{i))- Each 
successor, j g S{i) adjusts its phase according to its own 
edge specific PRC, fji{4>i) where — >■ + fji{(t>j) (si- 
multaneous signals are processed sequentially in random 
order). For example, in the well known model of Mirrolo 
and Strogatz [T] this phase adjustment rule is given by 
fij = y~^(e -I- V{<j)i)) — <j)i for concave V, while the ex- 
treme case where fij{4>i) = 1 — <j)i leads to a resetting of 
oscillator i with immediate firing while the other extreme 
case of fij{(f>i) — —(pi leads to a resetting of oscillator i 
without firing. The most well studied models of PCOs 
are either purely excitatory {fij{(pi) > 0) or purely in- 
hibitory (/„■((/),) <0). 

For the sake illustration, consider the PRCs we de- 
note "strong reseting" (SR), where for some Bp G (t, 1), 
ftji^'i) = ~<t>i for < 01 < So otherwise /^j = 0. Syn- 
chrony is clearly a solution for these curves; every oscil- 
lator is simply reset to time r after all oscillators fire. 
To study this solution consider the time 1 + t map, H . 

It turns out there is also a clear way to understand part 
of the basin of the SR system near synchrony. Denote 
(/i+(t) = maxi((/)i(t)), (t>^{t) = mini(0i(t)) and p{t) = 
(f>'^(t) — it). Furthermore let po{x, y) = min(a; — r, 1 — 
y -\- t), where x and y are some system parameter. Then 
in the SR case if at time t' no signals are en route and 
the range p{t') < po{Bo,Bq) then a careful analysis of 
the system shows that 

iJ(0j) = min(</ii + T,minjgp(i)((/)j))- 

Notice that if an oscillator j succeeds an oscillator with 
phases 4'~ then H{(j)j) — . In this way the minimum 
spreads, first to the successors of the minimum and then 
to the successors of the successors and so on. If the graph 
is aperiodic, then there exists some d such that d appli- 
cations of the successor function, denoted: S'^, includes 
the whole graph. Thus on aperiodic graphs this process 
leads to synchronization. 
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SR PRCs harness inhibition to stabihze synchrony, yet 
in many cases (such as when p > .5) excitation is use- 
ful, as it can aUow an oscillator to push forward a cas- 
cade. It is possible to augment SR PRCs with excitation 
while preserving similar convergence bounds and meth- 
ods of analysis, in the special case where the graph is 
undirected (directed aperiodic graphs are dealt with by 
a later theorem). Consider the PRC we denote "strong 
firing" (SF) where fij{(t>i) = ~(f>i for < cfii < Bq and 
otherwise fij{4>i) = 1 — 0i- Notice that in this situation 
oscillators always end up with phase after they receive 
a signal, either by being reset or if the phase is greater 
than Bq, by firing. 

To understand the map H in the SF case consider the 
excitation and inhibition separately, and again constrain 
p{t') < pq{Bo, Bq). It can be verified that between appli- 
cations of H every oscillator fires. Let Xi{to) be the next 
time that i fires after some time Iq. During an application 
of H excitation leads to: 

Xi{to) = min(to + 1 - ^^(io), minjgp(j)Aj(to) +t). 

Applying this map to an undirected graph implies that all 
oscillators fire within r of their predecessors. The effect 
of the inhibition then behaves very similar to that of the 
SR case, where the relative phase differences are captured 
by min(— Ai -|- t, minjgp(i) — A^). An interesting feature 
of the min map in the inhibition is that it preserves the 
T predecessor-successor phase difference, which implies 
that after a single iteration of H the oscillators will no 
longer receive signals in the excitatory regime, leading to 
the same type of convergence as with SR PRCs. 

These two phase response curves converge because of 
the way a modified min map spreads across a graph, di- 
rected aperiodic in the case of SR and undirected aperi- 
odic for SF. While the convergence of SR and SF PRCs 
is easy to understand other PRCs converge more ro- 
bustly (as will be demonstrated later) and require a 
more general argument. Consider the family of PRCs, 
"strong type 11" (S2) which have the following require- 
ments. The first requirement has an initial "reset zone" , 
fiji4'i) = for < < r -|- K, where k > 0. In 
addition the response curve must be slightly less than r- 
inhibitory, in that fij{4>i) < —t — k, on [r, Bq] and then 
be excitatory, fij{(j)i) > 0, for > Bi & [Bq,!). An 
illustration of S2 curves can be seen in Fig. [T] Notice 
that in {Bo,Bi) there are no restrictions on /y, though 
the smaller this region is the more general the conver- 
gence results will be. Similarly, it is not necessary that 
the PRC encode the value of r in its shape, only that it 
is greater than r inhibitory. 

We can now generalize the results for the SF and SR 
PRCs to this new class. Furthermore this result remains 
true even if the graph changes over time. Indeed, sup- 
pose that the graph changes every 1 + r time, giving 
graphs Gi, G2, . . . Gfc. This leads to our main result: If 
each Gfc has no isolated nodes, if there exists d such that 




FIG. 1. (Left) The value of r and the shape of S2 PRCs 
determine the values of -Bo and Bi and thus the strength of 
the convergence results. (Right) Displayed are three different 
empirically generated PRCs, where the red dashed curve is 
from ventricular heart cell in an embryonic chick [4] , the solid 
blue curve is from rabbit sinus node cells [5] and the data 
points with the line fit to them is from a low-threshold spiking 
GABAergic interneurons 

Si{Si+i{. . . Sd{v) ■ ■ ■)) — V for an infinite number of /, 
and if the initial range p < min(i?o — t, 1 — _Bi + t) 
the system will converge to synchrony. Furthermore if 
the previous conditions hold for all I then convergence 
occurs before time t* — pd/min(T, k). The graph condi- 
tions might at first seem onerous, but many such exam- 
ples abound. For example, a grid that suffers a random 
edge failure each time period would suffice, as would one 
where each Gfc is a different random tree. While the pre- 
vious results for SF and SR PRCs were based on an un- 
derstanding of how one phase spreads across the graph, 
this argument follows by focusing on the worst case sets 
of signals an individual oscillator can receive. 

First: if the initial conditions are contained in an in- 
terval of size p and there are no signals in transit, it is 
easy to see that this will remain true under iterations of 
the time 1 + t map. To see this, translate time so that 
the oscillators with the largest phase are just about to 
fire at time t — 0: 0~''(O) = 1^. The key insight is that 
all signals will occur within at most time p+r, since each 
oscillator can only fire once in that time, and any oscilla- 
tor that has not fired will be in the excitatory region. In 
addition, note that the oscillator with the largest phase 
can not receive a signal until at least time r so will al- 
ways be inhibited by at least r. Thus after time 1 -|- t it 
will be at most about to fire again, since it was inhibited 
at least r and never excited. A careful analysis of the 
remaining oscillators using these insights shows that the 
time 1 + T map does not increase the size of the interval 
of phases. 

Next, we consider oscillators with (j)i{Q) < 1 ~ e where 
e = min(K,T). Such an oscillator will not fire until at 
least time e. Now consider the successors of such an 
oscillator, j with (0) > 1 — e. Oscillator j will receive a 
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signal from ia,tt> e+r when it is in inhibitory region, so 
will be inhibited by a sufficient amount such that (f)j{l + 
r) < 1 — e. A careful argument also shows that if 4>j{0) < 
1 — e this will remain at time t+r. Iterating this argument 
d times (recall that 5''' is the complete graph) on the 
oscillator with the smallest phase, shows that the time 
+ t) map will reduce the size of the phase interval by 
at least e, or if it is less than e the phases will completely 
synchronize, proving our convergence result. 

The benefit of the more general class of phase response 
curves can be seen in Fig. [2] In Fig. [2]the PRC "Limited 
Resetting", which has values strictly greater than —.1, is 
more likely to converge than SR PRCs when run on a 
slightly modified binary tree with uniform initial con- 
ditions. Indeed the combination of inhibition and cuts 
that divide a graph into multiple disjoint subgraphs can 
lead to solutions where nodes along the cut never fire 
(and where the phases are never smaller than the critical 
p) . The risk of these non synchronous solutions increases 
with the amount of inhibition in the PRCs, the number 
of different disjoint subgraphs and the size of the cuts. 




strong Resetting 
Limited Resetting 



FIG. 2. Limited Resetting offers more robust convergence on 
a binary tree with a single triangle and uniform random initial 
conditions than SF. (The 'S' like shape of the graph is due to 
whether the system converges to a phase on the same side of 
the tree as the triangle or far.) 

An important application of this result arises in sensor 
networks [3, that is, collections of many small sen- 
sors which communicate over radio frequencies. There 
has been great interest in the use of PCOs to provide a 
simple and robust mechanism to synchronize sensor net- 
works. Such systems have intrinsic propagation and pro- 
cessing delays. In addition, most sensor networks have a 
complex graph structure, and a not necessarily constant 
network topology. However, most theoretical analyses ig- 
nore delays and assume the complete graph. As seen in 
Fig[3j these assumptions are quite significant. In partic- 
ular under the assumptions in our first result, the type 



II PRCs converge rapidly and robustly, while the cur- 
rently most popular PRCs [6] have bounded error but 
fail to synchronize exactly. In addition, even when we 
consider more realistic conditions, such as variations in 
propagation delay times and heterogeneous oscillator fre- 
quencies, the type II PRCs still perform well, while the 
top competitors do not. 

Furthermore, the system is relatively robust to error. 
For example. Notice that the provable basin of attraction 
is the greatest when i?o = = -5 + t, in which case the 
the system converges so long as p < .5. In this case, with 
probability 1 the synchronous solution is robust to any 
single random error in the oscillator phases. 

Our detailed convergence analysis allows us to "tune" 
the PRCs for specific objectives. For example, as dis- 
cussed in [5] in sensor-net applications it is important 
that the oscillators be allowed to "sleep" for as long as 
possible between firings. In our model this corresponds 
to choosing a large interval [Bq,Bi] and setting the re- 
sponse curve to in that region. However, choosing a 
large interval decreases the basin of attraction for the 
synchronous state. A simple and robust solution for this 
problem, would be to start with small interval [Bq^Bi] 
and then expand it over time. For example, one could 
start with Bq = Bi = 1/2 + t which allows for a large 
initial interval of phases and then reduce Bq (and increase 
Bi) by T or less every 2d{l + t) units of time, stopping 
when Bq = 2t. This provides provable convergence with 
a long sleep period outside of an initial period. 

In the case where the delay t is known in advance then 
the oscillators can simulate the arrival of their own signal 
to themselves r time after they fire. This modification 
has the effect of adding self loops to the graphs, vastly 
increasing the kind of underlying sensor networks that 
this system performs well on. Indeed, if ever node has a 
self loop the only additional requirement on the graph is 
that the graph never becomes permanently unconnected. 

Type II PRCs have also been seen in many places in 
nature |H] . As seen in Fig. [T] actual phase response 
curves taken from cells in the heart and from some cor- 
tical interneurons are described by S2 PRC suggesting 
that the stability of synchronous solutions in these set- 
tings may be described by our system. Furthermore, that 
this family of curves was discovered by a genetic algo- 
rithm [TO] lends credence to both the evolability and the 
performance of such PRC for providing synchrony. 

However, thus far the discussion has included only 
unweighted graphs, yet many biological graphs are 
weighted. An extension of this model that allows for 
weighted graphs is explored in [13j . This extension allows 
each edge's PRC to be weighted such that fij{(j)) = —(p 
for (j) < Wij and less than Wij for x G [wij , Bq] . In 
this situation convergence follows from the condition that 
^jWij > T. Thus generalizing the result to weighted 
graphs. 

This weighted version can also be used to give interest- 
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FIG. 3. (a) Under the theoretical setting with a uniform 
delay of 5% S2 curves outperform others from ChaosOS l6j, 
IEEE05 (6], and SIAM90 It (b) This trend continues under 
more realistic settings, with frequency error up to 2.5% and 
transmission and processing delay also up ot 2.5% (c) The 
random geometric graph on which these simulations was run. 

m 



gin and the only requirement on the magnitude of the 
inhibition is that the PRC have a slope of —1 at the 
origin and be nonzero in the inhibitory domain. 

One can also use high indegree to provide strong prob- 
abilistic convergence results. For example, for any S2 
PRC that is nonzero in the excitatory region there exists 
a constant c such that if the indegree of the connection 
graph is larger than clog(n/e) then under uniform ran- 
dom initial conditions the system will converge to syn- 
chrony with probability of at least e. In general, for a 
well chosen S2 phase response curve one can give explicit 
bounds on the probability of convergence for any non 
zero initial distribution of oscillator phases by analyzing 
the indegrees. 

Other modifications are possible too. One such mod- 
ification is the addition of "quiescent" period g, where 
oscillators ignore signals for time q after receiving a sig- 
nal. This modification allows for a more general class of 
PRCs that are allowed more freedom between and t. 
This new system can be shown to converge on all strongly 
connected graphs, even periodic graphs. Another mod- 
ification is to allow for the oscillator to adjust not only 
their PRC but also their frequency. Numerical results 
suggest that the proper choice of the frequency response 
curve allows for oscillators with a more robust region of 
convergence. 

In summary, the family of S2 phase response curves 
that was introduced is relatively general and includes 
curves that were empirically found in systems which syn- 
chronize in nature. Furthermore, we outlined a proof 
that this family of PCOs has a robust region in which 
it converges to synchrony on strongly connected aperi- 
odic graphs. This convergence remains even in the pres- 
ence of uniform time delay and particular mutations in 
the graph. It was then noted that in numerical trials 
this method outperforms similar PCO based time syn- 
chronization methods for wireless sensor networks. This 
advantage remained in numerical runs even when hetero- 
geneous time delays, frequencies, and random errors were 
introduced. 
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ing analytic results in the situation where we have some 
advance knowledge of the underlying graph topology. For 
example, if the indegree of the graph is known to be at 
least fc, then the PRCs only need to be "resetting" over 
the interval from to r//c. This implies that in the limit 
of high indegree, the resetting region approaches the ori- 
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